home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATHLIB2 / MATH387.PAS < prev    next >
Pascal/Delphi Source File  |  1995-10-14  |  11KB  |  303 lines
  1. Unit MATH387;
  2.  
  3. (* Bibliotheque mathematique pour utilisation du coprocesseur flottant
  4.   JD GAYRARD Sept. 95
  5.  la bibliotheque est batie à partir des fonctions du coprocesseur
  6.  du type 386, elle fournit les fonctions suivantes:
  7.       fsin, fcos, ftan, farctan, farctan2,
  8.       farcsin, farccos, fmod, mod_2PI,
  9.       ten_to, y_to_x, fexp, fln, flog...
  10.  
  11. Aucune verification du domaine de definition des fonctions n'est faite,
  12. pas plus qu'un controle de la validite des operandes. Il est conseille
  13. d'utiliser cette bibliotheque pour les types single et double exclusivement *)
  14.  
  15. {$G+}
  16. {$N+}
  17. {$E-}
  18.  
  19. { table opcode du 387 non comprise par turbo pascal V7 }
  20. { FSIN    : D9 FE
  21.   FCOS    : D9 FF
  22.   FSINCOS : D9 FB
  23.   FPREM1  : D9 F5 }
  24.  
  25. interface
  26.  
  27. const author  = 'GAYRARD J-D';
  28.       version = 'ver 0.0 - 10/95';
  29.  
  30. type float = double; { a modifier suivant l'utilisation }
  31.  
  32. (* use only with 86387, 86486 or pentium for type single, double andt extended,
  33. no check of definition domain of the function or range (FPU limitation).
  34. The f prefix avoids function redefinition of system runtime library *)
  35.  
  36. function fsin(x : float): float;
  37. function fcos(x : float): float;
  38. procedure sincos(x : float; var sinus, cosinus : float);
  39. function ftan(x : float): float;
  40. function farcsin(x : float): float;
  41. function farccos(x : float): float;
  42. function farctan(x : float): float;
  43. function farctan2(x, y : float): float;   { return arctan(y/x) }
  44. function fmod(x , y : float): float;      { return x mod y }
  45. function fmod_2PI(x : float): float;
  46. function fln(x : float): float;
  47. function flog(x : float): float;
  48. function fexp(x : float): float;
  49. function ften_to(x : float) : float;      { return 10**x }
  50. function fy_to_x(y , x : float): float;   { return y**x }
  51. function module(x, y : float): float;
  52.  
  53. procedure ssincos(x : float; var sinus, cosinus : single);
  54. procedure dsincos(x : float; var sinus, cosinus : double);
  55.  
  56. implementation
  57.  
  58. function fsin(x : float): float; assembler;
  59. {if x < pi.2^62, then C2 is set to 0 and ST = sin(x)
  60.                  else C2 is set to 1 and ST = x }
  61. {no check range validity is performed in this function}
  62. asm
  63.    FLD x           { load x }
  64.    DB $D9, $FE     { opcode for FSIN }
  65. end;
  66.  
  67. function fcos(x : float): float; assembler;
  68. { if x < pi.2^62, then C2 is set to 0 and ST = sin(x)
  69.                   else C2 is set to 1 and ST = x }
  70. {no range validity check is performed in this function}
  71. asm
  72.    FLD x         { load angle }
  73.    DB $D9, $FF   { opcode for FCOS }
  74. end;
  75.  
  76. procedure dsincos(x : float; var sinus, cosinus : double); assembler;
  77. (* retourne sinus et cosinus(x), utilisable uniquement
  78.  avec 80387, 80468 et pentium et type double *)
  79. asm                          { ST(0)     ST(1) }
  80.      FLD x                   {  x          -   }
  81.      DB $D9, $FB             { cos(x)   sin(x) }
  82.      LES DI,cosinus          {                 }
  83.      FSTP ES:QWORD PTR [DI]  { sin(x)      -   }
  84.      LES DI,sinus            {                 }
  85.      FSTP ES:QWORD PTR [DI]  {  -          -   }
  86. end;
  87.  
  88. procedure ssincos(x : float; var sinus, cosinus : single); assembler;
  89. (* retourne sinus et cosinus(x), utilisable uniquement
  90.  avec 80387, 80468 et pentium et type single *)
  91. asm                          { ST(0)     ST(1) }
  92.      FLD x                   {  x          -   }
  93.      DB $D9, $FB             { cos(x)   sin(x) }
  94.      LES DI,cosinus          {                 }
  95.      FSTP ES:DWORD PTR [DI]  { sin(x)      -   }
  96.      LES DI,sinus            {                 }
  97.      FSTP ES:DWORD PTR [DI]  {  -          -   }
  98. end;
  99.  
  100. procedure sincos(x : float; var sinus, cosinus : float);
  101. (* retourne sinus et cosinus(x), utilisable uniquement
  102.  avec 80387, 80468 et pentium *)
  103. var lcos, lsin : float;
  104. begin
  105.      asm                 { ST(0)     ST(1) }
  106.           FLD x          {  x          -   }
  107.           DB $D9, $FB    { cos(x)   sin(x) }
  108.           FSTP lcos      { sin(x)      -   }
  109.           FSTP lsin      {  -          -   }
  110.      end;
  111. cosinus := lcos;
  112. sinus := lsin
  113. end;
  114.  
  115. function ftan(x : float): float; assembler;
  116. { if x < pi.2^62, then C2 is set to 0 and ST = 1 and ST(1) = tan(x)
  117.                   else C2 is set to 1 and ST = x }
  118. {no range validity check is performed in this function}
  119. asm              { ST(0)    ST(1) }
  120.    FLD x         {  x        -    }
  121.    FPTAN         {  1      tan(x) }
  122.    FSTP ST(0)    { tan(x)    -    }
  123. end;
  124.  
  125. function farcsin(x : float): float; assembler;
  126. (* retourne l'arcsin de x *)
  127. {  methode :                ________
  128.    arcsin(x) = arctan( x / V 1 - x.x ) }
  129. {no range validity check is performed in this function |x| > 1 }
  130. asm                 { ST(0)     ST(1)     ST(2) }
  131.    FLD X            {  x         -         -    }
  132.    FLD ST(0)        {  x         x         -    }
  133.    FMUL ST(0), ST   { x.x        x         -    }
  134.    FLD1             {  1         x.x       x    }
  135.    FSUBRP ST(1), ST { 1 - x²      x        -    }
  136.    FSQRT            { sqrt(1-x²)  x        -    }
  137.    FPATAN           { arcsin(x)   -        -    }
  138. end;
  139.  
  140. function farccos(x : float): float; assembler;
  141. { retourne arccos(x)
  142.    methode :            ________
  143.    arcsin(x) = arctan( V 1 - x.x / x ) }
  144. { pas de controle de domaine de definition |x| > 1 }
  145. asm                 { ST(0)     ST(1)     ST(2) }
  146.    FLD X            {  x         -         -    }
  147.    FLD ST(0)        {  x         x         -    }
  148.    FMUL ST(0), ST   { x.x        x         -    }
  149.    FLD1             {  1         x.x       x    }
  150.    FSUBRP ST(1), ST { 1 - x²      x        -    }
  151.    FSQRT            { sqrt(1-x²)  x        -    }
  152.    FXCH             {   x         z        -    }
  153.    FPATAN           { arccos(x)   -        -    }
  154. end;
  155.  
  156. function farctan(x : float): float; assembler;
  157. asm              { ST(0)    ST(1) }
  158.    FLD x         {  x         -   }
  159.    FLD1          {  1         x   }
  160.    FPATAN        { atan(x/1)  -   }
  161. end;
  162.  
  163. function farctan2(x, y : float): float; assembler;
  164. { retoune arctan (y / x) }
  165. asm              { ST(0)    ST(1) }
  166.    FLD y         {  y         -   }
  167.    FLD x         {  x         y   }
  168.    FPATAN        { atan(y/x)  -   }
  169. end;
  170.  
  171. function fmod(x, y : float): float; assembler;
  172. (* retourne x mod y *)
  173. asm                 { ST(0)    ST(1) }
  174.    FLD Y            {  y        -    }
  175.    FLD X            {  x        y    }
  176. @repeat_mod:
  177.    FPREM            { x mod y   y    }
  178.    FSTSW AX
  179.    SAHF
  180.    JP @repeat_mod
  181.    FSTP ST(1)       { x mod y   -    }
  182. end;
  183.  
  184. function fmod_2PI( x : float): float; assembler;
  185. (* retourne x mod 2.pi *)
  186. asm                   { ST(0)    ST(1) }
  187.    FLDPI              {  pi        -   }
  188.    FADD ST, ST        { 2.pi       -   }
  189.    FLD x              {  x        2.pi }
  190. @unit_circle:
  191.    FPREM              { x mod 2pi  2pi }
  192.    FSTSW AX
  193.    SAHF
  194.    JP @unit_circle
  195.    FSTP ST(1)         { x mod 2pi   -  }
  196. end;
  197.  
  198. function fln(x : float): float; assembler;
  199. { retourne le logarithme naturel de x, utilise
  200.  la methode loge(x) = loge(2).log2(x)}
  201. { pas de verification du domaine de definition (x < 0) }
  202. asm             {  ST(0)          ST(1)  }
  203.    FLDLN2       { ln(2)            -     }
  204.    FLD X        {   x             ln(2)  }
  205.    FYL2X        { ln(2).log2(x)    -     }
  206. end;
  207.  
  208. function flog(x : float): float; assembler;
  209. { retourne le logarithme base 10 de x, utilise
  210.  la methode log10(x) = log10(2).log2(x) }
  211. { pas de verification du domaine de definition (x < 0) }
  212. asm             {  ST(0)          ST(1)  }
  213.    FLDLG2       { log10(2)         -     }
  214.    FLD X        {   x           log10(2) }
  215.    FYL2X        {log2(x).log10(2)  -     }
  216. end;
  217.  
  218. function fexp(x : float): float; assembler;
  219. { retourne e^x, par la methode e^x = 2^(y.log2(e))
  220. { 2^z = 2^f.2^i with f = frac(z) and i = int(z)
  221. { 2^f is computed with F2XM1, 2^i with FSCALE }
  222. const round_down : word = $177F;
  223. var control_ww : word;
  224. asm                   { ST(0)     ST(1)     ST(2) }
  225.    FLD X              {  x         -         -    }
  226.    FLDL2E             { log2(e)    x         -    }
  227.    FMULP ST(1), ST    { x.log2(e)  -         -    }
  228.    FSTCW control_ww
  229.    FLDCW round_down
  230.    FLD ST(0)          {  z         z          -   }
  231.    FRNDINT            { int(z)     z          -   }
  232.    FLDCW control_ww
  233.    FXCH               {  z         i          -   }
  234.    FSUB ST, ST(1)     {  f         i          -   }
  235.    F2XM1              { 2^f-1      i          -   }
  236.    FLD1               {  1        2^f-1       i   }
  237.    FADDP ST(1), ST    { 2^f        i          -   }
  238.    FSCALE             { 2^f.2^i    i          -   }
  239.    FSTP ST(1)         { 10^x       -          -   }
  240. end;
  241.  
  242. function ften_to(x : float): float; assembler;
  243. { retourne 10^x, par la methode 10^x = 2^(y.log2(10))
  244. { 2^z = 2^f.2^i with f = frac(z) and i = int(z)
  245. { 2^f is computed with F2XM1, 2^i with FSCALE }
  246. const round_down : word = $177F;
  247. var control_ww : word;
  248. asm                   { ST(0)     ST(1)     ST(2) }
  249.    FLD X              {  x         -         -    }
  250.    FLDL2T             { log2(10)   x         -    }
  251.    FMULP ST(1), ST    { x.log2(10) -         -    }
  252.    FSTCW control_ww
  253.    FLDCW round_down
  254.    FLD ST(0)          {  z         z          -   }
  255.    FRNDINT            { int(z)     z          -   }
  256.    FLDCW control_ww
  257.    FXCH               {  z         i          -   }
  258.    FSUB ST, ST(1)     {  f         i          -   }
  259.    F2XM1              { 2^f-1      i          -   }
  260.    FLD1               {  1        2^f-1       i   }
  261.    FADDP ST(1), ST    { 2^f        i          -   }
  262.    FSCALE             { 2^f.2^i    i          -   }
  263.    FSTP ST(1)         { 10^x       -          -   }
  264. end;
  265.  
  266. function fy_to_x(y, x : float): float; assembler;
  267. { retourne y^x, par la methode y^x = 2^(y.log2(y))
  268. {no range validity check is performed in this function (y > 0) }
  269. { 2^z = 2^f.2^i with f = frac(z) and i = int(z)
  270. { 2^f is computed with F2XM1, 2^i with FSCALE }
  271. const round_down : word = $177F;
  272. var   control_ww : word;
  273. asm                   { ST(0)     ST(1)     ST(2) }
  274.    FLD Y              {  y         -         -    }
  275.    FLD X              {  x         y         -    }
  276.    FYL2X              { x.log2(y)  -         -    }
  277.    FSTCW control_ww
  278.    FLDCW round_down
  279.    FLD ST(0)          {  z         z          -   }
  280.    FRNDINT            { int(z)     z          -   }
  281.    FLDCW control_ww
  282.    FXCH               {  z         i          -   }
  283.    FSUB ST, ST(1)     {  f         i          -   }
  284.    F2XM1              { 2^f-1      i          -   }
  285.    FLD1               {  1        2^f-1       i   }
  286.    FADDP ST(1), ST    { 2^f        i          -   }
  287.    FSCALE             { 2^f.2^i    i          -   }
  288.    FSTP ST(1)         { y^x        -          -   }
  289. end;
  290.  
  291. function module(x, y : float): float; assembler;
  292. (* retourne le module du complexe (x,y) *)
  293. asm                  { ST(0)     ST(1) }
  294.    FLD Y             {  y         -    }
  295.    FMUL ST(0), ST    { y.y        -    }
  296.    FLD X             {  x         y.y  }
  297.    FMUL ST(0), ST    { x.x        y.y  }
  298.    FADDP ST(1), ST   { d.d        -    }
  299.    FSQRT             {  d         -    }
  300. end;
  301.  
  302. end.
  303.